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Abstract 

A two-dimensional Navier-Stokes code has been 
developed for rapid numerical simulation of axisymmetric 
flow fields, including flow fields with an azimuthal veloc- 
ity component. The azimuthal-invariant Navier-Stokes 
equations in a cylindrical coordinate system are mapped to 
a general body-fitted coordinate system, with the stream- 
wise viscous terms then neglected by applying the thin- 
layer approximation. Turbulence effects are modeled 
using an algebraic model, typically the Baldwin-Lomax 
turbulence model, although a modified Cebeci-Smith 
model can also be used. The equations are discretized 
using central finite differences and solved using a multi- 
stage Runge-Kutta algorithm with a spatially- varying time 
step and implicit residual smoothing. 

Results are presented for calculations of supersonic 
flow over a waisted body-of-revolution, transonic flow 
through a normal shock wave in a straight circular duct of 
constant cross-sectional area, swirling supersonic (invis- 
cid) flow through a strong shock in a straight radial duct, 
and swirling subsonic flow in an annular-to-circular dif- 
fuser duct. Comparisons between computed and experi- 
mental results are in fair to good agreement, 
demonstrating that the viscous code can be a useful tool 
for practical engineering design and analysis work. 

Introduction 

Axisymmetric flow fields are frequently encountered 
in internal and external aerodynamics, within various 
engineering research areas such as ballistics, aircraft inlets 
and nozzles, and circular/annular internal duct flows. 
Often these flows involve velocities in the azimuthal 
direction. For the rapid simulation of such flow fields, a 
computer code has been developed which solves the thin- 
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layer formulation of the azimuthal-invariant Navier- 
Stokes equations, including the momentum equation in the 
azimuthal direction. 

The original work in developing and computationally 
applying an azimuthal-invariant thin-layer form of the 
Navier-Stokes equations, including the momentum equa- 
tion for the azimuthal direction, was done by Nietubicz et 
al. [1-2]. In that work, the three-dimensional, trans- 
formed, thin-layer equations in Cartesian coordinates were 
reduced to an azimuthal-invariant set by an appropriate 
coordinate selection and subsequent integration in the azi- 
muthal direction to determine axisymmetric source terms. 
Since the thin-layer approximation was applied before 
reducing the equations, only the in viscid source terms 
were determined. In the present work the thin-layer 
approximation has been applied after transformation of the 
axisymmetric Navier-Stokes equations to the computa- 
tional domain, so that viscous source terms were retained 
in the equation set. At moderate to high Reynolds num- 
bers, however, where the thin-layer approximation is 
valid, these viscous sources are probably always small and 
could be neglected. Nevertheless, they have not been dis- 
carded since their retention has little impact on computa- 
tional efficiency. 

The numerical solution of the governing equations is 
performed using an explicit time-marching scheme with 
central finite differences for spatial discretization. The 
time marching is carried out to steady-state conditions, 
with the convergence rate accelerated using a spatially- 
varying time step and implicit residual smoothing. The 
baseline code is robust, executes rapidly, and vectorizes 
well on vector processors. 

Several flow field solutions are presented that demon- 
strate the application, accuracy, and practical usefulness of 
the method. Specifically, results are presented for calcula- 
tions of supersonic flow over a waisted body-of-revolu- 
tion, transonic flow through a normal shock wave in a 
straight circular duct of constant cross-sectional area, 
swirling supersonic (inviscid) flow through a strong shock 
in a straight radial duct, and swirling subsonic flow in an 
annular-to-circular diffuser duct. Experimental or analyti- 
cal data are compared to the computed results, except in 
the last test case, with generally fair to good agreement 
between the two. 
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Governing Equations 

The Navier-Stokes equations are written in a station- 
ary cylindrical coordinate system, and all terms with 
derivatives in the azimuthal-direction are excluded. The 
cylindrical equations are mapped to a general body-fitted 
Q coordinate system using standard techniques. As 
illustrated in Figure 1, the ^-coordinate direction is along 
the body-surface, and therefore approximately along the 
flow direction, while the ^-direction is roughly normal to 
the flow. The viscous derivatives in the ^-direction are 
dropped by applying the thin-layer approximation [1-3]. 
Using ( x, y, z) for denoting the cylindrical coordinate 
directions, with x in the axial direction, y in the azimuthal 
direction, and z in the radial direction, the resulting equa- 
tions can be written as follows: 

d { q + d^ + d^[d-Re~ 1 ^) = fl + Re^K ( 1 ) 
where 
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The absolute velocity components u, v, and w point in the 
x, y, and z coordinate directions, respectively. The contra- 
variant velocity components are given by 

U=^ x u + % w 

(2) 

Assuming an ideal gas with constant specific heats, then y 
(which is c p /c v ) is constant, and the energy and static 

pressure are given by 

c = p [7<7n) + ^“ 2 * v2+ "' 2 )] (3) 

p = (r- 1) [« , -|p(« 2 + v 2 + >v 2 Jj (4) 

where the sonic velocity a is related to static pressure and 
density by the equation of state: 



The terms in the viscous flux vector S can be written 
as follows: 



where the stress terms are 

°xx = + 

+ W 0) 

a zz = + XV-V 

The terms in the viscous source vector K are given by 
*3 * 

( 8 ) 

k a = -XV- V 

Stokes' hypothesis, X = - 2p,/3, is assumed, and the 
dynamic viscosity |i and the thermal conductivity K are 
dimensionless. The velocity divergence V * V in the dila- 
tation term is given by 


v-v = + 



(9) 
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The inverse-Jacobian J 1 is given by the relation 

7 _1 = I = z(x^-x^) (10) 

with the metric terms defined as follows: 


inlet total conditions. The Reynolds number Re and 
Prandtl number Pr are defined in terms of these properties. 

For turbulent flows the dimensionless dynamic vis- 
cosity |ul is determined by 

» = ^lam + Vturb ( 16 ) 



(11) 


The metric terms and Jacobian are calculated from the 
grid coordinates by using second-order central 
differences. In order for the metrics to be geometrically 
consistent with the Viviand strong conservation form of 
the governing equations [4], the following constraints 
must be met: 


Vc = Vs 

These constraints can be satisfied in the discretized met- 
rics by using the following difference equations: 

^ >x ” 4 ( Z i t j+ 1 + Z iJ- 1^ ( Z i,j + 1 ” z ij- 1^ 

= l < z i,j + i +z .j-i> <• x ij + 1 - x ij - 1 } 

j (13) 

~^>x ~ 4 ^/+ l,j + Z i- \,P ^ z i + l,j~ 

^ = i (z i+u +z i-h? 

where the i- and y'-indices correspond to the and indirec- 
tions, respectively. In order to avoid numerical sources, 
the inverse-Jacobian is calculated using a rewritten form 
of Equation (10): 

J~ l = z[^[-ZJcJ + ^(^)] (14) 

with the radius z outside the brackets taken at the central 
node, and central differences applied to the discretized 
form of the terms in parenthesis: 


= ^ Z U + l +Z /,J-l) ( X i,j+ 1 ~ X i,j-\} 
= \ {z i+hj +z i-h? ( x i + i,r x i-iJ 


(15) 


The equations are nondimensionalized by constant 
reference properties; namely, density p r ^ , sonic velocity 

a re j , and dynamic viscosity , all typically based on 


where the laminar viscosity is calculated using a power- 
law function of dimensionless temperature T: 

M ta = (<r-I>!-)” = (a 2 )” (17) 

with n = 2/3 for air over a wide range of temperatures. 
The turbulent eddy- viscosity \^ turb is calculated using an 

algebraic turbulence model, either the Baldwin-Lomax 
model [5] or a modified version of the Cebeci-Smith 
model [6], The vorticity magnitude used in the turbulence 
model is calculated according to the following relation: 

M = Ju 2 x + (i) 2 y + m 2 z (18) 

where 

% = ^ + 

-m z = ^ + 

The dimensionless thermal conductivity term, K over 
Prandtl number Pr, is determined as follows for turbulent 
flows: 


JL 

Pr 


jf lam ^ turb 

P r lam P r turb 


(19) 


with the turbulent Prandtl number assumed equal to 0.90. 

The axisymmetric governing equations can be 
reduced to the planar two-dimensional set by eliminating 

the source terms H and K on the right-hand-side of Equa- 
tion (1), setting the azimuthal velocity v to zero every- 
where, setting the (undifferentiated) radius z equal to one 
in Equations (10) and (11), and eliminating the term w/z 
in Equation (9). 


Boundary Conditions 

Different boundary types can apply to a particular 
problem depending on the specified geometry and whether 
the problem involves an internal or external flow field. 
The boundary types addressed below include the axis-of- 
symmetry, the solid wall with and without the no-slip con- 
dition, and the inflow and outflow boundaries for internal 
flow fields. Boundary conditions suited mostly for exter- 
nal flow fields are not addressed. 

At the inflow boundary the total temperature, total 
pressure, and whirl (v-velocity) are fixed according to the 
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specified inflow condition. For flows with subsonic 
(throughflow) velocities at the inflow boundary, a steady 
characteristic relation using the upstream-running Rie- 
mann invariant is used to solve for the remaining quanti- 
ties at the boundary. The Riemann invariant is based on 
the meridional velocity component, and the meridional 
flow angles are fixed at the specified input values [7-9]. 
Since in internal flows the inflow boundary can involve 
large radial velocities and viscous shear layers at walls, a 
non-isentropic formulation of the steady, axisymmetric 
characteristic equation is applied. Assuming that the 
inflow direction is nearly parallel to the ^-direction, the 
steady axisymmetric characteristic equation can be written 
as follows: 



The Riemann invariant based on meridional velocity is 
given by 

R = V~— - (21) 

m y — 1 

where 

V = Ju 2 + w 2 (22) 

tn 

The flow variables on the right-hand-side of Equation (20) 
are taken from the interior, and all derivatives are obtained 
by first-order (backward) differences from boundary to 

interior. The resulting formulation is solved for R at the 
inlet boundary. Note that the entropy-gradient term T is 
non-zero only within viscous regions. The meridional 
velocity at the inlet is calculated from the total tempera- 
ture and the Riemann-invariant by using the following 
equation: 


where 


y _ (Y-l)/?~ + .yg-2C 
m 7 + 1 


(23) 
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4II+I) 
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(y+l)v 2 + (y-l) 



For convenience, the term B uses the dimensionless total 
sonic velocity rather than the total temperature. Once the 
meridional velocity is known, the velocity components u 
and w can be algebraically determined using the specified 
meridional flow angle. Finally, using isentropic flow rela- 


tions the density is calculated from the total velocity mag- 
nitude and the specified total temperature and total 
pressure. 

For supersonic inflow all flow quantities at the inlet 
boundary are fixed. 

At the outflow boundary the conservation variables p, 
p w, pv, and pw are extrapolated from the interior. For flow 
fields with supersonic outflow the energy conservation 
variable e is extrapolated as well. For flow fields with 
subsonic outflow the static pressure p is specified at the 
outer radial position, and the remaining local static pres- 
sures along the exit boundary are calculated by integrating 
the simple radial-equilibrium equation: 

•*P = (24) 

dz z 

Along solid walls with a slip-condition (for inviscid 
computations) the flow is obtained by linear extrapolation 
of the density p, the static pressure /?, the azimuthal veloc- 
ity v, and the covariant velocity U ' given by 


IT = 


X^U + 

r 2 2 




(25) 


The wall tangency condition is satisfied by setting the con- 
travariant velocity W equal to zero. 

The axis-of-symmetry boundary condition is similar 
to the slip-wall condition except that the radial velocity w 
is identically zero. The azimuthal velocity v is also zero 
for viscous computations, while for inviscid computations 
it is linearly extrapolated. The axial velocity u and static 
temperature T are determined using second-order one- 
sided differences such that = 0 . The static 

pressure is linearly extrapolated, and the density follows 
from the equation of state. Note that for inviscid computa- 
tions the axis flow quantities are determined only for aes- 
thetic reasons since they do not influence the flow field 
solution. 

For solid walls with the no-slip condition (for viscous 
computations) the velocities u and w are zero, and 
v = £lz , where £2 is the angular velocity of the wall. The 
wall static pressure is obtained by zero-order extrapolation 
from the adjacent interior point; that is, = 0 . The wall 

static temperature is either specified on input, or else 
determined assuming an adiabatic wall, 7^ = 0, and 

using a second-order one-sided difference. 


Multistage Runge-Kutta Scheme 

The governing equations are discretized using a node- 
centered finite-difference scheme, with second-order cen- 
tral differences used throughout. The multistage Runge- 
Kutta scheme developed by Jameson, Schmidt, and Turkel 
[10] is used to advance the flow equations in time from an 


4 



initial guess to a steady state. Equation (1) can be rewrit- 
ten in the following terse form: 

d t q = -J[R r (R v +D)} (26) 

where R f is the inviscid residual including source term, 

is the viscous residual including source term, and D is 

the artificial dissipation term described in the next 
subsection. This allows the multistage Runge-Kutta algo- 
rithm to be expressed as follows: 

*0 = % 

— ^-VArtfl^o ~(Ry + D)q (>] 

(27) 

9% = <7o” vtfcJAt 1 — <*v + ^ id 

4n+ 1 = % 

In most applications a standard four-stage scheme is used, 
with 


a * = [l/4, 1/3, 1/2, l] 

For efficiency both the physical and artificial dissipation 
terms are calculated only at the first stage and then held 
constant for subsequent stages. 


Artificial Dissipation 

The dissipative term D in Equation (26) is similar to 
that used by Jameson et al. [10]. It is given by 

Dq = {D^ + D^q (28) 

where the ^-direction operator is given by 

D % q = ^ [T* C 5 ( V 2 \q + V^q) ] (29) 

The term is a coefficient, and the terms V 2 and V 4 are 
the following: 


V 2 = \i 2 msa{v i _ v \ i ,v i + v v i+2 ) 
V 4 = max (0, \i A -V 2 ) 

where 


(30) 


|Pi + 1 -2p, + p,_il 

^{p n \p i + l + ^P i +P l . l \) 


(31) 


and the i-indexing corresponds to the ^-direction. The 
constant | x 2 scales a first-order artificial viscosity that is 

switched on at shocks detected by Equation (31). The 
denominator in Equation (31) is normally constant at the 


pressure p n , making the operator roughly symmetric 

across shocks. For external flows p n is equal to the static 

pressure at infinity, and for internal flows it is normally a 
fixed low pressure at either the grid inlet or exit. The more 

commonly applied term | P i+ \+ 2p. + p,-_ j| is included 

to switch on the second-difference dissipation when the 
pressure becomes very small, usually due to numerical 
problems. The constant p 4 scales a uniform third-order 

artificial viscosity that is switched off at shocks by Equa- 
tion (31). Typically, values for \i 2 and 1 1 4 are 1/8 and 

1/64, respectively. 

The coefficient C ^ is arbitrary, but can have a large 

impact on stability and solution accuracy. It has been 
found in general that directionally-biased coefficients 
work substantially better than directionally-independent 
ones. A directionally-biased coefficient which works well 
is similar to one proposed by Kunz and Lakshminarayana 
[11], which was a modification of an earlier one proposed 
by Martinelli and Jameson [12]. For the two-dimensional 
axisymmetric case it can be written as follows for each 
grid direction: 


* 

C c = ( 


At r 

1 + ^ 




A u 

1 + ^ 


(32) 


. —1 . -1 
A \ + 


To minimize artificial dissipation in viscous regions, 
the coefficients in the ^-direction are reduced to zero qua- 
dratically over several grid points near walls. The one- 
dimensional time steps in Equations (32) are given below 
in connection with the stability limit. 


Stability Limit 

The following expression, obtained in part through a 
linear stability analysis, and in part by transforming three- 
dimensional Cartesian linear stability results [8,9,13] to an 
axisymmetric coordinate system, is used for the time step: 

* 

— zi — it < 33 > 

+A +A t a +A t v 

* 

where X is the maximum Courant number for the partic- 
ular multistage scheme, which for the standard four-stage 
scheme is approximately 2.8. The inverse one-dimen- 
sional time steps are given by 
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< 



. -1 
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+ a j^X + ^Z 


II 

1-^ 
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II 


viscous time-step contribution At 1 , the c 


(34) 


is given a value of 4.0, and an axisymmetric term xy is 
included where 



To accelerate convergence to a steady state, the maxi- 
mum permissible time step at each grid point is used giv- 
ing a constant Courant number everywhere. The time step 
is normally updated every 10 iterations. 


Implicit Residual Smoothing 

To further accelerate convergence it is desirable to use 
a time step even larger than the stability limit given by 
Equation (33). To maintain stability, the residual calcu- 
lated in Equation (26) is smoothed after each Runge-Kutta 
stage by an implicit smoothing operator: 

(l-e^)(l = R k (36) 

where 5^ and 6^ are standard second difference opera- 
tors, and and are smoothing parameters. 

Linear stability analysis shows that the Runge-Kutta 
scheme may be made unconditionally stable using implicit 
residual smoothing if the smoothing parameters 8 are 
made sufficiently large [14]. In one dimension 



gives unconditional stability if X is the Courant limit of 
the unsmoothed scheme, and X is a larger operating Cou- 
rant number. In two dimensions different e’s may be used 
in each direction, and their magnitudes may often be 
reduced below the value given by Equation (37). Good 
success has been achieved by using the artificial dissipa- 
tion coefficients and in Equations (32) to scale e in 

the two directions: 




C^E 

C c e 


(38) 


which is an approach similar to that of Martinelli and 
Jameson [12]. Typically, the value for e can be taken as 
the minimum limit from Equation (37), which usually 
works well for Courant numbers X up to 7.0. 


Computational Details 

The code used to simulate the flow fields discussed 
below is referred to as DVC2D (Duct Viscous Code 2D), 
and it was written primarily for internal flow fields. The 
code can be run as an inviscid (Euler) solver, or as a vis- 
cous solver with or without turbulence modeling. In all 
viscous cases presented below, the Baldwin-Lomax turbu- 
lence model was used and updated every five iterations. 

Computational grids were generated using an interac- 
tive program called TIGGERC [15] in conjunction with 
interior grid restretching using an elliptic partial differen- 
tial equation (PDE) solver [16]. The elliptic PDE solver 
forces grid orthogonality at solid boundaries, and allows 
user control of the grid spacing normal and adjacent to the 
boundaries, and of the stretching away from the bound- 
aries. 

The results discussed in Reference 6, along with some 
recent computational studies done by the authors, reveal 
that an adjustment to the Baldwin-Lomax turbulence 
model parameters C Kleb and C Qp gives better solutions 

for turbulent flat-plate boundary layers. The suggested 
parameter values are 


C Kleb “ 0 646 

C = max (1.216, 0.8 Af ) 

cp v e ' 


(39) 


where M e is the local Mach number at the edge of the 

boundary layer. Notice in Equation (39) that C Qp has a 

constant value of 1.216 for local Mach numbers less than 
1.52. The adjustment for higher Mach numbers is only a 
rough, first-order approximation and was determined from 
computations of turbulent boundary layers on a flat plate 
with freestream Mach number 2.4. Comparisons between 
experimental [17] and computed velocity profiles in the 
wake region of the boundary layer indicated that a value of 
1.920 was significantly better, at free-stream Mach 2.4, 
than the 1.216 recommended for subsonic (incompress- 
ible) boundary layers [6]. The comparison is shown in 
Figure 2 where the distance from the wall y is nondimen- 
sionalized by the momentum thickness 0. The general 

effect of increasing C cp at the higher Mach numbers is to 

increase the turbulent eddy viscosity in the outer region, 
resulting in a thicker boundary layer and larger skin fric- 
tion. 

Normally 2000 to 3000 time-step iterations, involving 
3 to 5 orders of magnitude reduction in the density residu- 
als, are needed to achieve convergence to a steady-state 
viscous solution. The execution time of the DVC2D code 
on different computer systems was measured for a super- 
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Figure 2 — Turbulent flat-plate boundary-layer velocity 
profile for freestream Mach 2.4 (experimental 
data reproduced from Reference 17) 


sonic viscous flow simulation on a 150 x 75 grid. Based 
on 500 iterations, the execution times for different com- 
puters systems are listed in Table 1. 

Core memory requirements for the code are small by 
current standards, being less than one megaword of ran- 
dom access memory (RAM) for a typical viscous simula- 
tion. 


Results 

Waisted Body-of-Revolution 

The first set of results to be presented involves calcu- 
lations of supersonic flow over a waisted body-of- 
re volution. The experiment is described in detail in Refer- 
ence 18, so only a brief discussion introducing the basic 
features of the experiment will be provided here. The 
waisted body-of-revolution is shown in Figure 3, and as 
described in the reference was designed to provide certain 
desired variations in the streamtube area and streamwise 


Table 1 — Approximate code execution times for various 
computer systems 


Computer system 
(CPU) 

time/node/iteration 

(ns) 

CRAY C-90 

3.6 

CRAYY-MP 

7.5 

SGI a Power Challenge 
(R8000 — 75MHz IP21) 

10.6 

SGI 4D/440 [4 CPU’s] 
(R3000 — 40MHz IP7) 

61.4 

SGI lndigo2 

(R4000 — 100MHz IP22) 

83.1 

SGI 4D/35 

(R3000 — 36MHz IP12) 

176. 


a Silicon Graphics, Inc. Compiler options were selected for high 
levels of optimization and varied with the workstation model. 


pressure-gradient, for the purpose of studying their influ- 
ence on turbulent boundary-layer character and develop- 
ment. 

Nominal test conditions at zero angle-of-attack 
included five different supersonic upstream Mach numbers 
ranging between Mach 1.4 and 2.8, as well as two sub- 
sonic Mach number conditions, Mach 0.6 and 0.8. The 
Reynolds number was varied for some tests, but measure- 
ments at the nominal freestream Reynolds number Re L of 

l.OxlO 7 based on the body length of 1.524 m (5 ft.) were 
made at each Mach number condition. Measurements of 
wall static pressure and local skin-friction were obtained 
along the length of the body, and boundary layer profiles 
were measured at a few selected axial locations. 

Only two of the test conditions were examined com- 
putationally for the present work; namely, nominal 
upstream Mach numbers 1.4 and 2.8, with freestream Rey- 
nolds number Re L equal to l.OxlO 7 . The same computa- 
tional grid was used for simulations at both conditions, 
and it is shown in Figure 4. The grid dimensions are 
150 x 75, but in the figure only every other inline (vertical 
line) is shown. The grid spacing normal and adjacent to 
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Figure 3 — Waisted body-of-revolution 
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Figure 4 — Computational grid for waisted body-of-revo- 
lution (only every other vertical grid line 
shown) 


Contour increment 0.05 



Figure 5 — Mach contours for computed flow field with 
upstream Mach 1 .398 


ized by body length, which is needed to imbed 3 to 4 inte- 
rior points within the viscous sublayer. 

The numerical simulation for upstream Mach number 
1 .4 was run assuming a turbulent boundary layer along the 
entire body, which is consistent with the experimental 
intent. As reported in Reference 18, a boundary layer 
transition trip of size 0.13 mm was attached 38.1 mm from 
the nose of die body. This was adequate for the Mach 1.4 
condition, but the trip proved to be too small for the higher 
Mach number conditions as discussed in the reference, and 
later here with respect to the Mach 2,8 results. 

Mach contours for the computed nominal Mach 1.4 
case are shown in Figure 5, where it is indicated that the 
upstream Mach number is, more precisely, 1.398. Only 
about half of the computational domain is shown since the 
outer region is not important for the comparisons to be 
made. The outer (far-field) boundary was treated as an 
inviscid reflecting surface, but its distance from the body 



Location x/L 


Figure 6 — Measure and computed pressure coefficient 
distributions for upstream Mach 1 .398 



Figure 7 — Measure and computed skin-friction coeffi- 
cient distributions for upstream Mach 1 .398 


is far enough away to avoid interference from the reflected 
waves. 

The measured and computed distributions of pressure 

coefficient C and of skin-friction coefficient C, are com- 
P J 

pared in Figures 6 and 7, respectively. As can be seen the 

agreement is good, especially for pressure coefficient. 

Note that the skin-friction coefficient is based on local 
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dynamic pressure, rather than upstream dynamic pressure 
as might be expected [18]. 

Measured and computed boundary-layer profiles in 
terms of Pitot pressure are compared in Figure 8. In 
Reference 18 the profile data are given in terms of velocity 
ratios and Mach numbers, but to eliminate the added 
uncertainties inherent in these quantities due to data reduc- 
tion assumptions, the measured Pitot pressures were 
recovered by reverse-calculation from the reported data. 
The agreement is generally good, although in the adverse 
pressure-gradient region, at x/L = 0.550 and 0.700, sig- 
nificant differences exist, as might be expected consider- 
ing the well-known limitations of the turbulence model. 

Mach number contours for the computed nominal 
Mach 2.8 case are shown in Figure 9, where as indicated 
the upstream Mach number is actually 2.799. Again, only 
about half of the computational domain is shown since the 
outer region is not important for the comparisons. The 
computed pressure coefficient distribution is compared 
with the measured distribution in Figure 10, and can be 
seen to agree well with it. 

The measured and computed Pitot pressure profiles 
for the boundary layer at several axial locations are shown 
in Figure 11. In this case two computed solutions are 
shown; the dashed lines corresponding to an all-turbulent 
boundary layer, and the solid lines corresponding to a 1am- 
inar/turbulent boundary layer (and to the computed results 
in Figures 9 and 10) with transition occurring at the onset 
of the adverse pressure gradient, at about x/L = 0.42 . 
As can be seen, the computer simulation with boundary- 
layer transition is in much better agreement with the data, 
confirming that the transition trip was too small at this 
Mach and Reynolds number. As discussed in Reference 


Contour increment 0. 10 



Figure 9 — Mach contours for computed flow field with 
upstream Mach 2.799 



Location x/L 

Figure 10 — Measure and computed pressure coefficient 
distributions for upstream Mach 2.799 
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Figure 11 — Measured and computed Pitot pressure profiles at several axial locations, for upstream Mach 2.799 


18, this was established after the initial experiments, and a 
trip size on the order of 0.33 mm would have been needed 
to force transition at the nose of the body. Furthermore, 
because the boundary layer was not fully turbulent the 
skin-friction data are of questionable accuracy, at least 
over the front half of the body — the measurement tech- 
nique was calibrated for only turbulent flow. For this rea- 
son, skin-friction results were not compared in this case. 

Closer examination and comparison of the measured 
and computed profiles in Figure 11, particularly the pro- 
files at x/L = 0.400 , reveals that the front half of the 
boundary layer probably was not completely laminar in 
the experiment. Some indications of turbulence are appar- 
ent in the shape of the experimental profile, but the turbu- 
lence may have been from the ffee-stream, rather than that 
arising from a transitional layer. In fact, a measurable 
influence of ffee-stream turbulence would be expected for 
a laminar boundary-layer in such a strong favorable pres- 
sure gradient. A more definitive evaluation cannot be 
given, however, since ffee-stream turbulence levels for the 
wind tunnel were not reported with the experimental 
results. 

Finally, some information pertaining to both the Mach 
1.4 and 2.8 numerical simulations might be included. 
First, as shown in Figure 12 the solutions behaved well 
numerically, with a fairly continuous decrease in the den- 
sity-residuals by over four orders of magnitude in 2000 
iterations. Both cases were run with a CFL number of 7.0, 
with residual smoothing coefficients at minimum values, 
and with low artificial viscosity levels — the first- and 
third-order artificial dissipation coefficients p 2 (scaled by 

1/4) and p 4 (scaled by 1/16) were equal to 0.50 and 


0.20, respectively, for the Mach 1.4 flow, and 0.50 and 
0.30, respectively, for the Mach 2.8 flow. The fully turbu- 
lent cases were simulated by setting the turbulence model 
transition parameter C mutm [5] equal to zero, whereas the 

Mach 2.8 flow with transition was simulated by setting 
C mu tm ec l ua l t0 20 — the normal (default) value of 14 

produced a “transitional” layer with the turbulence model 



Figure 12 — Convergence histories for waisted body-of- 
revolution flow field simulations 


10 




turning briefly on, then off through the favorable pressure- 
gradient region, and then on again as the adverse pressure- 
gradient was encountered. 

Straight Circular Duct 

The next set of results to be presented involves the 
calculation of transonic flow through a straight circular 
duct of constant cross-sectional area. The experiment is 
described in detail in Reference 19 and was intended for 
studying the flow field of a normal-shock-wave/turbulent- 
boundary-layer interaction. It was performed in a test sec- 
tion consisting of a tube with inside-diameter D of 0.2477 
m, and length L of 2.690 m, and with the inside surface 
fabricated so as to be hydraulically smooth. Supersonic 
flow with a nominal Mach number of 1.5 entered the tube 
from a supersonic nozzle, and a normal shock wave was 
positioned within the tube at various locations. A hollow 
cylinder inserted into the exit of the tube was used to 
establish and position the shock wave. Because of wall 
boundary-layer growth a slight Mach-number gradient 
existed along the length of the test section, leading to some 
variation in the shock-incident Mach number as the shock 
position changed. 

The test condition considered for the present work 
was one at which the most extensive measurement series 
was carried out [19-21]. It had a reference Mach number 
of 1.44, at Reynolds number Re D of 5.43xl0 6 based on 

the inside diameter of the tube — the reference point is at 
the boundary-layer edge just upstream of the shock wave. 
At this test condition both unsteady and mean-flow (slow- 
response, steady-state) measurements were performed in 
the test section, although only the mean-flow data are used 
here. Mean-flow measurements included wall static pres- 
sures, and probe-surveys with Pitot pressure, static pres- 
sure, and total-temperature probes. Mean-flow shear- 
stresses at the wall were determined from measurements 
with a Preston tube. 

It is important to note that the survey and wall shear- 
stress data were obtained by positioning the shock wave at 
various locations in the test section, relative to a single, 
fixed measurement location. The wall static-pressure dis- 
tribution relative to the shock wave could be obtained in 
the same way; that is, by using a single static-tap at the 
probe location. However, it could also be obtained for a 
fixed shock position, or even several fixed shock positions, 
by using the array of static pressure taps distributed along 
the length of the test section. The different possibilities 
this allows seem to account for variations in the wall 
static-pressure distributions presented in the literature [19- 
21] over a period of several years, for the same experi- 
ment. 

In view of the experimental approach, comparisons 
between measured and computed results are somewhat 
problematic. There is a significant quantitative difference 
between the (shock-relative) flow field for a fixed shock 
location, and that for a variable shock location. Neverthe- 
less, comparisons are made, but with care in attempting to 
draw valid inferences. 


The numerical simulation of the flow field was sensi- 
tive to small changes in flow conditions. This was 
expected for a transonic flow field with strong viscous/ 
inviscid interaction, where even small changes in 
upstream or downstream conditions would substantially 
alter the shock position. Consequently, it was important 
that the supersonic upstream boundary condition to the 
computational domain closely match that of the 
experiment. This could be done for only one of the exper- 
imental shock positions, which was the position where the 
shock wave was located farthest downstream. 

A few assumptions and simplifications were involved 
in using the measured shock-incident flow profile to deter- 
mine the upstream boundary condition for the computa- 
tional domain. Briefly, the turbulent boundary layer 
profile was constructed using the measured skin-friction 
coefficient of 0.0018 in conjunction with the mea- 
sured velocity profile data [21]. Near the wall, for dimen- 
sionless wall distances y + less than 50, Spalding’s inner- 
law composite formula [22] was used to define the veloc- 
ity profile. Between that region and the innermost mea- 
sured velocity data points, a cubic-spline interpolation of 

dimensionless velocity m + versus In (^y + J was used to 

fill in the rest of the profile. The static pressure was 
assumed constant at the reference value across the bound- 
ary layer, and the temperature profile was approximated by 
assuming an adiabatic wall and using the Crocco-Buse- 
mann [23] approximation: 


T 

T 


T 2 00 




v 


(40) 


where T aw is the adiabatic wall temperature, and the tur- 
bulent recovery factor r was given a value of 0.895. The 
resulting density profile agrees well with the measured 
density profile, as shown in Figure 13, where the distance 
from tiie wall y is nondimensionalized by the tube inside- 
diameter D. (The resulting velocity profile will be shown 
later along with several downstream profiles.) The profile 
in the free-stream, between the boundary-layer edge and 
the centerline, was constructed using the measured veloc- 
ity profile in conjunction with the measured plenum total 
temperature and the reference static pressure, assuming 
homentropic flow. 

The computer simulation was performed on two dif- 
ferent grids, both of dimensions 150 x 50, with the same 
grid spacing from wall-to-centerline. The grid spacing 
normal and adjacent to the wall is about 2.0x1 0“ 5 , as non- 
dimensionalized by tube diameter, which imbeds at least 3 
to 4 interior points within the viscous sublayer. The 
upstream halves of the two grids are shown in Figure 14 
along with Mach contours for the corresponding computed 
solutions. The so-called “coarse grid” has uniform spac- 
ing in the axial direction, whereas the “fine grid” has been 
locally (and manually) refined around the normal shock 
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Figure 13 — Measured and specified density profiles at 
inlet boundary with reference Mach 1 .44 


Figure 15 — Measured and computed pressure coefficient 
distributions for reference Mach 1 .44 
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Location (x-x Q )/D 

Figure 14 — Computational grids and Mach number con- 
tours for computed flow fields 


for better shock resolution. The coarse grid is considered 
more typical of what would be used in routine code appli- 
cation where the location of strong shock waves is often 
not fixed, nor known a priori. 

To obtain the solutions in Figure 14 it was necessary 
to iteratively vary the exit static pressure until the shock 
was located just downstream of the inlet boundary. How- 
ever, the shock was not pushed so far upstream that the 
beginning of the shock/boundary-layer interaction region 
(at x = x Q ) reached the inlet boundary. In both cases the 

shocks are close enough to the inlet boundary that the flow 
profile does not change significantly between inlet and 
shock, as can be inferred from the Mach contours. 

The simulations were run with a CFL number of 7.0, 
with residual smoothing coefficients at minimum values, 
and with (scaled) first- and third-order artificial viscosity 
coefficients of 0.50 and 0.25, respectively. Due to the 
strong viscous/inviscid interaction and marginal steady- 
state flow-field stability, around 20,000 iterations were 
required to achieve full convergence (5 to 6 orders of mag- 
nitude reduction in the density-residuals) at a particular 
exit pressure. Such a large number of iterations is very 
abnormal and should be considered peculiar to this type of 
problem. 

Comparisons between the measured and computed 
wall static-pressure distributions are shown in Figure 15. 
The two sets of experimental wall pressure data corre- 
spond to two different measurement approaches; one with 
the shock wave at a single location [20] in the test section 
(diamonds), and the other with the shock wave relocated 
to various positions [21] relative to the probe-survey loca- 
tion (circles). The case with a single shock position is, in 
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principle, most like the computation, but for that case the 
exact inlet flow profile is not known (and therefore not 
used for the computation). In general, the computed dis- 
tributions for the coarse and fine grid are almost identical 
to each other, and are in basic agreement with the experi- 
mental data. 

The measured and computed skin-friction coefficient 
distributions are compared in Figure 16. As with the wall 
static-pressure results the two computations are almost 



Figure 16 — Measured and computed skin-friction coeffi- 
cient distributions for reference Mach 1 .44 


identical to each other, but here their agreement with 
experiment is not as good. The agreement might, how- 
ever, be considered fair from the standpoint of accuracy 
requirements for much engineering work, and in view of 
the strong adverse pressure-gradient. Results very similar 
to this were reported by Viegas and Horstman [20], and 
the failure of the computation to more accurately simulate 
the flow field lies with the turbulence model. Reasons for 
the lack of accuracy are well-known and understood in the 
field of turbulence-modeling research and development; 
for example, see Reference 24. 

Notice that the computed solutions indicate a rela- 
tively small, weak separation “beneath” the shock wave, 
that is, between ( x - x Q ) /D equal to 0.2 and 0.6, with a 

rapid increase in skin friction once reattachment occurs. 
The experimental results seem inconclusive in this regard, 
with the original study indicating that there was separation 
[19], but later work suggesting that separation probably 
did not occur [20]. 

Velocity profiles at several axial locations are com- 
pared in Figure 17. Also shown in the figure is the 
upstream profile, at x = x Q , which is a fit to the experi- 
mental data as explained earlier. Consistent with the static 
pressure and skin-friction results just presented, the agree- 
ment is generally fair between the measured and computed 
velocity profiles. 

A basic summary of the results helps to clarify the 
comparison. First, the computed flow in the interaction 
region exhibits an initial compression (pressure rise) 
which is too rapid, and which is followed by too gradual of 
a subsonic flow diffusion (pressure rise). The differing 
rates of subsonic diffusion can be seen in the wall static- 
pressure gradients, as well as in the progression of velocity 



Figure 17 — Measured and computed velocity profiles at several axial locations, for reference Mach 1.44 
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profiles (see Figure 17, where the experimental free- 
stream velocity level decreases more rapidly with distance 
downstream than does the computed velocity level). The 
measured values of skin-friction are higher than those 
computed, consistent with and corresponding to a larger 
decrease in boundary-layer-wake momentum, i.e., lower 
velocities. 

In Figure 17, one particular feature is worth further 
discussion, namely, the condition at the boundary-layer 
edge for the profile just downstream of the shock wave, at 
(x-x^)/D = 0.41 . At this location the computed 

velocity (fine grid solution) at the boundary-layer edge is 
slightly supersonic at about Mach 1.02 ( u/a { = 0.96), 

compared to a measured value of Mach 1.04. This com- 
puted (see Figure 14) and measured region of supersonic 
flow at the edge of the boundary layer just downstream of 
the normal shock, sometimes called a “supersonic 
tongue”, is mentioned in Reference 19, and has been 
observed by others for similar flows [25-27]. The super- 
sonic tongue results from the bifurcated, or “lambda”, 
shock structure at the base of the normal shock, and appar- 
ently exists only with shock-incident Mach numbers 
greater than about 1.4. Within the lambda shock system 
the compression occurs more gradually, as a supersonic 
diffusion, with less total-pressure loss. 

Straight Radial Duct 

The next simulation case to be presented involves 
swirling supersonic flow through a strong shock in a 
straight radial duct. The primary purpose of this test case 
is to demonstrate the application of the method to an axi- 
symmetric flow field with swirl, while at the same time 
allowing an assessment of the solution accuracy under 
such conditions. Documented experimental data for tests 
involving compressible swirling axisymmetric flow fields 
are difficult to find, so an inviscid “two-dimensional” sim- 
ulation has been run for which there is an exact analytical 
solution to compare with. A set of solutions involving a 
strong shock in the duct was selected mostly for the sake 
of added interest, but also to demonstrate the accurate 
shock-capturing capability in the radial direction. 

The geometry and grid for the simulation can be sim- 
ply described and so they need not be shown. Briefly, the 
duct is straight and radial, extending from a dimensionless 
radius of 1.0 to a radius of 2.0, with an aspect ratio 
(A z) / ( Ajc) of 10. The grid had dimensions 150 x 15, 
evenly spaced in both directions (15 nodes in the ^-direc- 
tion, or ^-direction, from wall-to-wall). The simulated 
flow field was two-dimensional with a zero velocity com- 
ponent in the jc-direction; that is, the inflow was a swirling 
radial flow with no variation wall-to-wall. 

Four different solutions, corresponding to three differ- 
ent subsonic exit static pressures and one all-supersonic 
case, were computed. The flow at the grid inlet was fixed 
at a radial Mach number component of 1.05, with azi- 
muthal (swirl) component 0.95, giving a resultant inlet 


Mach number of 1.416 at a direction 42.14 degrees from 
radial. The simulations were run for 2000 iterations with a 
CFL number of 7.0, with (scaled) first- and third-order 
artificial viscosity coefficients of 1.1 (up to 1.5) and 0.50, 
respectively. The residual smoothing coefficients in the £- 
direction were at minimum values, while some additional 
smoothing in the ^-direction was necessary to obtain good 
convergence behavior — a continuous decrease in the den- 
sity-residuals by about five orders of magnitude. 

The computed results are compared with the exact 
solutions in Figure 18, where graphs of Mach number, 
flow angle, and total pressure ratio are shown. Except in 
close proximity to the shocks, the agreement is virtually 
“exact”, as might be expected. Even through the shock 
waves, however, the agreement is very good in terms of 
both shock location and shock strength. Oscillations up- 
and down-stream of the shocks are minimal for a numeri- 
cal scheme employing central-differences, although to 
avoid stronger oscillations it was necessary to increase the 
first-order artificial viscosity coefficient somewhat (from 
1.1 to 1.5) as the shock strength increased. 

Annular-to-Circular Diffuser Duct 

The last simulation case to be presented involves 
swirling subsonic flow in an annular-to-circular diffuser 
duct. In this case no comparison is made with experiment 
or theory, the main purpose being to demonstrate the 
method for a swirling flow with viscous effects. Although 
it would be preferable to have a comparison with experi- 
mental data, as already mentioned it is difficult to find doc- 
umented data for experiments involving compressible 
swirling axisymmetric flows. 

The geometry with grid for the annular-to-circular 
duct is shown in Figure 19, where locations of three axial 
stations downstream of the centerbody are also indicated. 
The outer radius z Q of the duct is constant, the inner-outer 
wall radius ratio at the duct inlet is 0.50, and the center- 
body is a cylindrical section of length x/z 0 equal to 0.50 
followed by a tangent ogive with length-to-diameter 
aspect ratio 3.80. The aft end of the ogive is spherically 
rounded with an end-radius r/z G of 0.066. The dimen- 
sions of the grid are 1 15 x 65, with a spacing normal and 
adjacent to the walls (and centerline) of about 5.0xl0“ 5 , as 
nori-dimensionalized by outer radius, imbedding 2 to 3 
interior points within the viscous sublayer. 

For the simulation, the flow at the duct inlet was a 
swirling free-vortex with a mid-height Mach number of 
0.405 and a flow direction 28.7 degrees from axial. Tur- 
bulent boundary-layer profiles were specified for the inner 
and outer walls of the inlet, each layer with a thickness of 
20 percent of annulus height. The downstream static pres- 
sure at the outer wall was 93.5 percent of the inlet total 
pressure, and based on the mid-height inlet Mach number 
and duct outer-diameter, the Reynolds number for the sim- 
ulation was 4.44x1 0 6 . 
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Figure 18 — Exact and computed solutions for swirling 
inviscid flow through a strong shock wave in 
a straight radial duct 



Figure 19 — Geometry and computational gird for annu- 
lar-to-circular diffuser duct 


Contour increment 0. 05 





Figure 20 — Contours of various flow quantities for the 
computed flow field 


Contour plots for the simulated flow field are shown 
in Figure 20 , and graphs of axial and azimuthal velocity 
versus radius are shown in Figures 21 and 22 , respectively, 
for the three axial stations indicated in Figure 19 . A dis- 
cussion summarizing these results follows. 

To begin, it might be pointed out that except for the 
wall boundary-layer regions and the region behind the 
base of the centerbody, the flow field remains closely irro- 
tational (homentropic free- vortex). Some diffusion occurs 
in the outer half of the duct, as expected due to the area 
increase, whereas in the inner half there is no diffusion. 
The lack of a diffusion along the duct inner half is due to 
the decreasing radius which causes the azimuthal velocity 
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Axial Velocity u/a tQ 


Figure 21 — Radial distributions of axial velocity at axial 
stations downstream of centerbody 



Azimuthal Velocity v/a fQ 


Figure 22 — Radial distributions of azimuthal velocity at 
axial stations downstream of centerbody 


to increase considerably, azimuthal velocity being 
inversely-proportional to radius. That is, the flow pro- 
gressing downstream toward the base of the centerbody 
remains part of a ffee-vortex field — notice the horizontal 
orientation of the azimuthal- velocity contour lines, in con- 
junction with an almost uniform total-pressure field — 
except within the wall boundary layer. 

Far more intriguing than the above, however, is the 
flow at and behind the base of the centerbody. At the aft 
end of the centerbody an abrupt change in the flow can be 
seen to occur, downstream of which a large total-pressure 
deficit exists, especially near the centerline. This transi- 
tion is ostensibly a form of vortex breakdown, the occur- 
rence of which is to be expected under the flow conditions 
simulated. In discussing various types of vortex-break- 
down transitions, Keller [28] has written the following 
with respect to this particular type: “It can be observed in 
isolation by directly generating a supercritical hollow-core 
vortex in an annular vortex tube with a suddenly ending 
center body.” Vortex break-down phenomena and theory 
are well beyond the scope and intent of this work, but a 
few descriptive comments related to the simulated flow 
solution might be of general interest. 

References dealing with the subject of vortex break- 
down [28-31] describe it as a wave-like phenomenon, 
analogous to abrupt-transition phenomena in compressible 
gas flows and in open-channel (gravity, free-surface) 
flows. In the numerical test case under consideration, the 
vortex-breakdown transition occurs between the so-called 
supercritical vortex flow state upstream of the centerbody 
base, and a subcritical vortex flow state downstream of the 
centerbody. Thus, it is analogous to a strong shock wave 


in a compressible gas flow, or to a hydraulic jump in a 
shallow open-channel water flow. The result of the transi- 
tion is, locally, a substantial decrease in Mach number, azi- 
muthal velocity, and total pressure (see Figure 20). 
Although not shown, a sudden jump to large amounts of 
vorticity also accompanies the transition, with the azi- 
muthal velocity distribution near the centerline approach- 
ing that of a solid-body rotation. The profiles in Figure 22 
help to show this, as well as the reduction in azimuthal 
velocity levels below what would exist for an irrotational 
flow field. 

The axial velocity field downstream of the vortex- 
breakdown transition is difficult to discern from the above 
figures only, so in Figure 23 a velocity-vector plot is also 
provided. As can be seen in Figure 23, or even partially in 
Figure 21, the stream wise flow stagnates in a “bubble” just 
behind the transition. The velocities within the bubble are 



Figure 23 — Velocity vector field in the vortex-breakdown 
region 
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still relatively large, however, but only in the azimuthal 
direction (see Figure 22, station 1). In effect, a portion of 
the flow just swirls around the centerline without going 
anywhere, while producing a large area-blockage to the 
axial flow. 

To conclude this last case, a few comments about the 
numerical simulation can be added. The simulation was 
run with a CFL number of 5.0, with residual smoothing 
coefficients somewhat above minimum values. The 
(scaled) first- and third-order artificial viscosity coeffi- 
cients had a value of 0.50. Relatively low meridional 
Mach number levels for the flow field (weak compressibil- 
ity) caused slower-than-normal numerical convergence, 
especially near and on the centerline where the flow field 
developed somewhat slower. The solution achieved full 
convergence (4 to 5 orders of magnitude reduction in the 
density -residuals) at around 4,500 iterations. The simula- 
tion was somewhat temperamental regarding the selection 
of numerical parameters, apparently due to the “viscous” 
grid spacing at the centerline, and perhaps aggravated 
slightly by the unusually high vorticity concentrated near 
there. As a result, some additional residual smoothing was 
necessary for stability. An additional point might be men- 
tioned regarding this; namely, that stability problems were 
encountered in some subsonic simulations where initial 
azimuthal velocities were not near quasi-steady equilib- 
rium values. Under such conditions, even when a con- 
verged solution was achieved the number of iterations was 
abnormally large. A convenient way to systematically cir- 
cumvent this problem was found, and involves initially 
setting the azimuthal velocities everywhere to zero except 
at the upstream boundary. As the flow field subsequently 
develops the circumferential momentum convects down- 
stream without causing numerical problems. 

Concluding Remarks 

A two-dimensional Navier-Stokes code has been 
developed for rapid numerical simulation of axisymmetric 
flow fields, including flow fields with an azimuthal veloc- 
ity component. The azimuthal-invariant Navier-Stokes 
equations in a cylindrical coordinate system are mapped to 
a general body-fitted coordinate system, with the stream- 
wise viscous terms then neglected by applying the thin- 
layer approximation. Turbulence effects are modeled 
using an algebraic model, typically the Baldwin-Lomax 
turbulence model, although a modified Cebeci-Smith 
model can also be used. The equations are discretized 
using central finite differences and solved using a multi- 
stage Runge-Kutta algorithm with a spatially- varying time 
step and implicit residual smoothing. 

Results were presented for calculations of supersonic 
flow over a waisted body-of-revolution, transonic flow 
through a normal shock wave in a straight circular duct of 
constant cross-sectional area, swirling supersonic (invis- 
cid) flow through a strong shock in a straight radial duct, 
and swirling subsonic flow in an annular-to-circular dif- 
fuser duct. Comparisons between computed and experi- 
mental results were in fair to good agreement, 


demonstrating that the viscous code can be a useful tool 
for practical engineering design and analysis work. 

The code might also find useful application as a plat- 
form for turbulence model development and assessment. 
Along this line, some near-term future development of the 
code is planned. Specifically, recent efforts by Chima [32] 
in implementing the k-co two-equation turbulence model 
[33] will be incorporated into the DVC2D code. The 
increase in execution time for the code is estimated to be 
around 15 percent, and core-memory requirements will be 
extended by about 30 percent. 

As a further enhancement, some investigation is cur- 
rently underway to examine the feasibility of incorporat- 
ing a preconditioning scheme into the solver [34-36]. This 
feature would improve the code considerably by extending 
its applicability to flow fields with very weak 
compressibility. A current limitation of the code, common 
to many codes of this type, is that convergence problems 
become acute when relatively large regions in the compu- 
tational domain are nearly incompressible. The presence 
of high azimuthal velocities has been observed to provide 
no relief in this regard; that is, compressibility as intro- 
duced solely through circumferential momentum is appar- 
ently insignificant with respect to the influence of 
compressibility on numerical convergence. 
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